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Report 

The activities of this grant are divided into two tasks, each with a graduate 
student, a faculty advisor, and a separate topic. 

For the first task, Mr. Mehtab Pervaiz is working with Professor Judson R. 
Baron on spatial and temporal adaptive algorithms for reacting flows. A pa- 
per was prepared for and presented at the First World Congress on Compu- 
tational Mechanics held at University of Texas at Austin in September 1986. 
The attached copy of the paper summarizes the research in progress. Between 
September and February, interest has shifted to two dimensional flows. 

For the second task, Mr. Steven Ruffin is working with Professor Earll M. 
Murman on hypersonic viscous calculations for delta wing geometries. The 
conical Navier-Stokes equations model has been selected in order to investigate 
the effects of viscous-inviscid interactions. The more complete three-dimensional 
model is beyond the available computing resources. The flux vector splitting 
method with van Leer’s MUSCL differencing is being used. Both two and three 
stage Runge-Kutta time marching schemes have been investigated. A detailed 
von Neumann stability analysis has been performed for selecting the coefficients 
for each stage. Preliminary results have been computed for several conditions. 
Mr. Ruffin visited NASA Langley in January 1987 to work for a week with Dr. 
James Thomas on development of the numerical approach. 
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1 Introduction 

This paper discusses the numerical integration of quasi-one-dimensional 
unsteady flow problems which involve finite rate chemistry and are expressed 
in terms of conservative form Euler and species conservation equations. 

The coupled behavior between fluid flow and finite rate chemistry can 
introduce appreciable stiffness into numerical integration schemes, which 
then involve prohibitively long computation times. The use of globally fine 
grid resolution to ensure the capture of local flow features can also result in 
lengthy runs. The aim of the present paper is to provide a description of a 
controlled grid resolution approach in both space and time, and to demon- 
strate its advantage in diminishing the stiffness constraint. 

The use of uniform spatial and temporal grids demands some form of 
equilibrium limit modeling for extremely fast kinetics. However, the reten- 
tion of a fine grid resolution when approaching that limit is required only 
for small portions of the overall space/time domain. Typically, fine spatial 
resolution is desired in those regions in which a shock or species concentra- 
tion results in very rapid local changes; similarly high temporal resolution 
is needed where there are large non-equilibrium source terms which produce 
large temporal gradients. We discuss an adaptive technique which refines the 
spatial and/or temporal grid whenever preselected gradients exceed certain 
threshold levels. In general the resolved grid field is itself unsteady, and its 
rate of change may vary from large values for certain unsteady problems to 
zero for steady state solutions. For example, the adapted grid for a moving 
shock must track the discontinuity. 

The present algorithm involves periodic examination of the evolving solu- 
tion, detection of those regions in which large spatial non-uniformities exceed 
a threshold limit, and subsequent subdivision of the corresponding grids. 
Reverse embedding (collapse) to a coarser mesh is allowed up to the initial 
(coarsest) global grid. Consistent pre-embedding in accord with the initial 
flow field is appropriate. 

For unsteady flow the temporal gradients must be monitored so as to 
maintain sufficiently small time-steps for adequate local resolution and sta- 
bility. However, temporal adaptation does allow for spatial variation of the 
cell time-steps, without which a global minimum time-step would apply and 
could be extremely costly. 
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2 Governing Equations 


The unsteady, conservation form, quasi-one-dimensional Euler equations, 
for a reacting gas may be represented by the vector equation 


Here 
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where p, u, e, p and A are the density, velocity, total energy per unit volume, 
pressure and stream-tube area, respectively. Y, denotes the mass fraction of 
the s tH species (« = 1, • • • , S) and W, is the associated chemical source term. 
The vectors U, F and W each have 3 + 5 components, however £f =1 Y, = 1 
so that only 5 — 1 species equations need to be considered. An equation of 
state for an ideal gas mixture is assumed 


p^pRT^Y./m. (3) 

$ 

Here m, is the molecular mass of species s and Z is the universal gas constant. 
The total energy is 


;=e ua ,. + r c - d 7 " i + t « 

P , ■'To. 2 

where Hj , is the heat of formation at reference temperature 7o e and C v . is 
the constant volume specific heat, assumed to be constant for the adaptive 
examples to be shown. 

For R reactions of the form 

s s 

X>, r A. ^ E&rA. , r = 1, • J2 (5) 

f=l 1=1 

stiffness arises from the source terms 

R 

w. = Am, £(/?„ - «.r) [n /r - n„] (6) 
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in accord with differences in the flj for forward, backward progress rates. For 
the rth reaction 

n, = K, n (*£)' P) 

,Vi 

where a = a, r ,P„ for / = /r, br. The rate constants are assumed to be of 
the generalized Arrhenius form 

Ki = AiT 1 exp{-E t /T) (8) 


3 Integration Scheme 

The integration basis is a generalization of the second order accurate 
Lax-Wendroff finite volume scheme in Ref. [1]. Variables U, fluxes F and 
the source terms W are stored at the nodes. The generalization introduces 
non-uniform spatial grids and chemical source terms. Consider the cells B 
and C adjacent to node j (Fig. 1). The temporal change in state at node j 


ItJ 

6Uj = U” +1 - U” = 

f " At+ ^[ A< ’ +0 < A,S) 

(9) 

or, using Eqs. (1), 



6U, = {W- d -f x W + ^ 


(10) 

where, for example, 

c 

II 

£ 

(11) 

The flux balance for cell C, e.g., yields the cell change 


A U c = W' c At c + & - F k ) 

(12) 

in which W£ may be modelled as am average 


W^iWj + W,) /2 

(13) 


Alternatively, for a more accurate contribution to node j use can be made 
of a A Ujc based on choosing Wq = Wj t in which case 


A U jC = W, At c + (Fj - F k ) — ^ (14) 

t±xc 
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We shall use Eq. (12) for the determination of At and Eq. (14) for the 
integration scheme. 

In terms of a non-uniform grid parameter 


_ A xb — A xc 
; Axb + A xc 

a second order accurate expression for the rate of change of a scalar variable 
<p is 



d# 

dx 


= sS < «* ■ - « *> + 0 < Al » Al *> 


(16) 


Using this and Eq. (14) to discretize Eq. (10) the resulting overall change 
6Uj is comprised of contributions from cells B and C, 


SUj = SUjB + 6Ujc 


(it) 


where 


IV, a = ^ [AO,* + (AF„ + £-AF,„ + *,) + *»A W jB ] 
IV jc - '-J 1 [a Ujc - ^ (aFo - jjJjA Fjc + * c ) + ^A W iC ] 
Here, for example, 

dF I 


(18) 


dF 

A F c = — 

du 


A Uc , A Fjc = 




AUyc 


(19) 


An artificial viscosity contributes to the above two changes and is of the form 


*c = (°j + <’k)(U j -U k )/4 

The nodal viscosity coefficients are obtained from 

Oj = Omin + S(k B + K C ) 

where, for example, 


Kc = 


Pj-Pk 


( 20 ) 


( 21 ) 


( 22 ) 


\Pj + Pk\ 

The constant 6 is chosen so that Oj 6 ]» typically between 0.05 

and 0.5. 
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4 Spatial Adaptation 

For a system of K equations each and every state variable U(k),k = 

1 ,•••,#, need not be considered to define the necessary spatial and tem- 
poral resolution. Pre-selected variables for such space and time resolution 
criterion are referred to here as spaerit and temperit , respectively. A typical 
choice for spaerit is density since it appears as a factor in each element of 
the state vector; furthermore, density differences are present for most flow 
fields including contact discontinuities, and have been shown to be the an 
appropriate choice for a class of problems by numerical experimentation, 
Ref. [2]. 

When carrying out computations a distinction must be made between 
interests in steady and unsteady situations. Steady state problems may 
involve multiple grids and other acceleration techniques, whereas unsteady 
flows may not involve such procedures. Spatial adaptation can be done 
occasionally in steady flows and the number of such operations equals those of 
the spatially embedded levels desired; however for unsteady flows the spatial 
adaptation may have to be applied after each time-stride. The latter is 
defined as the time interval between the constant time surfaces for the entire 
spatiotemporal domain. It is more appropriate to use local time-stepping 
than temporal adaptation for steady flows. 

The choice of initial grid conditions is of some importance for unsteady 
flow. If large spatial gradients are present in the initial flow field and the 
spatial grid is coarse in their vicinity, the initial integrated solution will 
be degraded and will propagate as such to other spatial locations at later 
time levels. For that reason all initial flow fields are examined for first dif- 
ferences of spaerit, and wherever the critical limit is exceeded the grid is 
pre-embedded, i.e., some spatial resolution is added at appropriate locations 
before integrating the equations. For example a shock tube problem suggests 
that finer cells be inserted near the contact discontinuity surface. 

Upon initiating the integration scheme, the evolving solution is period- 
ically reexamined for regions of relatively large gradients of spaerit; if the 
threshold limit is exceeded the grid is divided locally . Alternatively, when 
associated gradients diminish on a previously refined grid, and become less 
than another critical limit, those contiguous grids may be collapsed while 
making certain that the cells to be merged are those from the same parent 
cell. The initial (coarse) global grid is kept stagnant by insisting that the 
coarsest cells (spatial level zero) be never merged to a coarser state, no mat- 


5 


ter how smooth the evolving solution proves to be. The choice of spatial 
threshold limits is similar to that for perfect gases [2,3]. 

The spatial interfaces are defined to be the nodes which are at the bound- 
ary of widely disparate cell sizes (|e| > |). A special formulation of artificial 
viscosity is used at the spatial interface for the steady flow situations due 
to the coupling with multiple grids. When integrating the coarse cells adja- 
cent to the spatial interface, the artificial viscosity coefficient is obtained by 
using Eq. (21) in the usual manner; however, for the neighbouring fine cell 
this coefficient is doubled when applying the contribution to the co mm on 
node for steady flow situations. Such a special formulation is unnecessary 
for unsteady cases since multiple grids are not used. 

Currently there are two major approaches for spatial adaptation proce- 
dures: grid point redistribution (or moving elements) and locally embedded 
grids. The current paper belongs to the latter class. A study that is typical 
among numereous grid point redistribution schemes is reference [4] while a 
local spatial embedding example is [5]. 

5 Temporal Adaptation 

The temporal and spatial adaptation procedures are inherently different 
and can be applied separately. The spatial adaptation process is carried out 
at the current time level at all spatial locations and may not be done as 
many number of times as the temporal adaptation. The frequency of spatial 
adaptation depends upon the time rate of change of the flow feature being 
resolved. The temporal adaptation, on the other hand, is repeated after 
each time-stride at all spatial locations and must anticipate the subsequent 
changes in the the flow field. The three steps for completing a temporal grid 
adjustment are : 

1. a determination of the allowable At for each cell 

2. reassignment of consistent At distributions such that adjacent cells 
differ by a factor of 2 at most 

3. determination of a proper integration sequence over the cell domain. 

The first order change with time from Eq. (12) is given for cell C by 

/ p F \ dU 

AU C = Ate (w; + -^) „ —A,C (23) 
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in term® of the west, east fluxes F w , F e ; and is subject to the usual CFL 
restriction 

At, t* < 

o/ being the local frozen speed of sound. In addition, the resolution re- 
quirement, A t rM , follows from a threshold criterion which is based on the 
maximum allowable change, AUc m „, for the tempcrit variable; i.e. 

\AUc\ < AI Jc mMl = \eiU c + *o| (25) 

where the e, sure small positive numbers. Effectively, the change is limited 
to a fraction of the state value excepting for vanishingly small levels. Unlike 
the choice of spaerit that for tempcrit is not obvious here. It should, how- 
ever, correspond to the maximum rate of change of a species. The fact that 
exchange reactions are generally faster than the dissociation reactions cam 
be important in selecting tempcrit. Comparing Eqs. (23) and (25) 

_ jci Uc + colA xc 
rM W c Ax c + (F„ - F e ) 

The resolution requirement, A t rei , may or may not exceed the stability re- 
quirement, A t # rt, and the actual time-step is 



A x c 
|u| + a, 


(24) 


At = min{At r „, At, t k} (27) 

Note that for non-reacting uniform flow At re , — ► oo and the stability re- 
quirement is governing. On the other hand, large Wc implies A t re , Atm, 
and the expected problem of stiffness. If time accurate descriptions are 
sought, a criterion such as Eq. (26) is required for both explicit and implicit 
integration schemes. Equation (26) also satisfies the stability restriction im- 
posed on the explicit schemes due to the chemical source terms. 

Based on the global minimum and maximum At constraints over the 
spatial domain, a current allowable temporal level m for the cells is assigned 
such that 

2 m < < 2 m+1 (28) 

Here m is equal to or less than a pre-selected maximum level, M, and is cho- 
sen to avoid what will be referred to as temporal level stiffness. A flow feature 
generally implies an associated characteristic speed, e.g. that for a shock or 
reaction, and the spatially embedded region must be sufficiently large to 
ensure that the feature will remain within the embedded region during the 
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subsequent time-stride. The size of a time-stride depends upon the value 
of M. An excessively luge size makes possible a departure of the feature 
from the spatially embedded region. Hence the embedded portion must be 
enlarged for large time-strides. Though temporal embedding involving large 
time-strides helps to reduce CPU time, the increased number of nodes can 
be expensive and a balance between the two competing effects is necessary. 
This represents one facet of temporal level stiffness. 

The actual time-steps are re-assigned for a given cell C, according to 

At% w = 2" A U (29) 

where level n < m is given by 

2" < m in{-^-,2 m } < 2 n+1 (30) 

The total number of time-strides for level n cells is 2 m-n . Figure 2 shows an 
example of cell time-steps before and after reassignment. The numbers in 
part (a) indicate the initial time-steps in terms of At min] part (b) shows the 
adjusted time-steps for m = 2. The temporal interface will be referred to as 
nodit (Node Of Different Time-steps) in the following. 

The order in which the integration over the cells takes place is of special 
importance. For example, if the level 0 cells labelled A< in Fig. 2b are 
integrated four times consecutively use is made only of information based on 
the two nodes d and e. This is correct for Ai, but for A4 additional account 
must somehow be taken of the nearby nodes. If for seven integration passes, 
we integrate level 0 cells on pass 1, level 1 cells on pass 2, level 0 cells on 
pass 3, level 2 cells on pass 4 and so on as indicated by the numbers in Fig. 2b, 
by the time A 4 will be integrated the nodes d and e will have accumulated 
effects from nodes a through <7, provided that after each integration pass 
the cells at a particular level have been updated and the flux, source terms 
and Jacobians recomputed. This represents a second facet of temporal level 
stiffness. 

The total number of integration passes is 

P T = 2 m+1 - 1 (31) 

A complete cycle of Pt integration passes will be referred to as a time-stride , 
and m its size. On pass P € (l, Pt] we integrate level n cells if 

p - 2" 

Jn^ T = integer (32) 
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No special formulation is needed at the nodits (which are not necessarily 
spatial interfaces) and in order to render the actual spatial location of any 
temporal level cell irrelevant, the data base must be constructed so as to store 
the cells with same temporal level together. There is no such restriction for 
the spatial adaptation pointers. 

Since both nodes of a cell are updated after each integration pass, and 
the flux, etc., are recomputed, the state at a nodit during a time-stride 
is not time-accurate; however, on completion of a time-stride the state is 
temporally correct. As an example consider the situation depicted in Fig. 3 
where Ate = 2A t$. On the first pass over cell B, the contribution to node 
j, SUjB, is computed from Eq. (18) based on the flux values at time level 
(n),so that a new state (*) at node j can be defined as 

u; = t/; + 6U*£ (33) 

The supercript (n, n) indicates that both nodes i and j of cell B use values 
at time level (n). The state (n + 1) at node * is defined in the usual manner 

U? +1 = U? + 6U& n + 6U& n (34) 

On the second pass cell C is integrated using values at the time level (*) for 
node j and time level (n) for node k and a new state (f) is defined at node 
j as 

uj = V) + SU'ff (35) 

where, SUJq is evaluated from Eq. (18), making use of, for example, 

A u;s = W7Alc + ( F i'~K)^ (36) 

Finally, on the third pass cell B is integrated again using values at time level 
(n + 1) for node i and values at level (f) for node j. The result is 

U? +2 = Uj + 6U% xa (37) 


6 Results 

Three examples have been considered to illustrate the technique: a stream- 
tube flow and shock tube, each with a single dissociating gas, and a multi- 
ple reaction flow for a diverging channel. For the first of these a Lighthill 
ideal dissociating gas [6], Mj *=* 2Af, was assumed to be flowing through 
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a converging-diverging streamtube with an area distribution A = 1 + 0.5z*. 
The nonequilibrium source term may be written as 

(38) 

with the non-dimensional reaction parameter given by 

♦ = C,V >' L ’ (39) 

m M yJp r /pr 

where Cf, tj , 6d and po are constants and subscript r denotes reference 
quantities. The rate parameter $ varies from zero for frozen flow, to infinity 
for equilibrium flow. 

An initial validation of the code consisted of shock free flow examples and 
comparison of results with [7] for several values of $. Figures 4 and 5 show 
results obtained with local time-stepping and a uniform grid; specifically the 
degree of dissociation and temperature distributions appear on a plot folded 
about the minimum area section such that the upper curves correspond to 
the subsonic upstream region. The freezing phenomenon is evident in the 
supersonic regions for intermediate $. The agreement with [7] was very 
good. 

Figure 6 shows a real flow through the same parabolic nozzle for $ = 10 4 
but with a back pressure ratio p»/p r = 0.92, so that a normal shock would 
be stationed at z = 0.5 for a frozen flow situation. Two levels of spatial 
embedding and local time-stepping were used for the adapted case. The 
results are shown corresponding to coarse, embedded and fine grids, with 
relative computing times 10.1 (fine/coarse) and 1.4 (adapted/coarse); the 
vertical scale corresponds to the coarse grid and the other two curves are 
displaced by the indicated offset. The embedded and fine grid solutions 
agree very well whereas the shock location is displaced and spread out for 
the coarse grid. Shown in Fig. 7 is the final grid and degree of dissocation 
for adapted and coarse grid cases, and indicates that the coarse grid solution 
predicts a different chemistry aft of the normal shock. 

A second example was the shock tube flow for frozen and reacting cases, 
($ = 0, 10 4 ). The initial conditions across the contact surface were pt/pi = 
0.2, Tj/Ti = 1.0. Three levels of both spatial and temporal adaptations 
were introduced and the results shown correspond to t = 0.6. Figure 8 
shows the progression of grids as time increases for the frozen case, whereas 
Fig. 9 shows the evolving temporal grid near time levels, t = 0 and t = 0.2. 


w M = *r»*> 
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Table 1 : Comparison of CPU time for shock tube 



$ = 0 

10 4 

10 ® 

Coarse 

1.00 

1.21 

1.89 

Adapted 

6.54 

8.77 

14.40 

Fine 

49.95 

65.38 

110.81 


Figure 10 indicates the computational results for $ = 10 4 at the final time for 
coarse, adapted and fine grids and Fig. 11 indicates the evolution of density 
on the adapted grid. We again note that the coarse grid solutions are poorer 
than the fine or adapted grid solutions. Although results corresponding to 
$ = 0 , 10 6 are not shown here, their CPU time comparisons are indicated in 
Table 1 . 

As a final example, multiple reactions were considered for the model in 
[ 8 ]. The reactions correspond to hydrogen combustion in air and are given 
by 

H 2 + 0 2 

Hi + 2 OH 

The forward rate constants (Eq. 8 ) are 


20H 
2 HjO 


(40) 


Ax 
A i 


(sijss _ 8<917 0 -28.95 ) x 10 44 m* / (kmolc.s) 
(l^s — o.833<£ -2.00 ) x 10 6 ® m 6 / {kmole 2 .s) 


TJl — — 10 T)i — — 13 

Ei = 2448.4 K E 2 = 18940.6 K 


where <j> is the equivalence ratio. Backward rate constants were formed by 
using the equilibrium constants [9] 


K ei = 26.164e~ 8992/r 

K t , = 2.682 x i 0 - 9 re +69415/r m 3 /kmole 
The streamtube area distribution is 



( 41 ) 
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For this case use was made of two levels of spatial and ten levels of temporal 
embedding. The calculations were carried out to steady state (rms error 
of mass-fraction of Hj < 10“* ) and took 4259s on a Microvax-II. For a 
source implicit scheme in [9] the same flow took 2524s on a CYBER-175 to 
converge. The comparative fine grid solution took three orders of magnitude 
longer. The results are shown in Figures 12 through 14 and proved to be in 
very good agreement with those of [9j. 

7 Discussion 

A strategy has been developed for automatic spatial and temporal grid 
embedding for a quasi-one-dimensional, inviscid reacting flow. The spatial 
resolution was added prior to the execution of each temporal cycle (time- 
stride), was based on first differences of the density field, and avoided spatial 
grid stiffness by limiting the non-uniformity parameter to 2:1 grid ratios 
for any adjacent cells. The procedure allowed for both refinement of and 
return to coarser grids, within some specified coarsest global spatial grid. 
The temporal adaptation allowed for a factor of two or more difference in 
local time-steps for contiguous cells. The larger the global disparity of the 
cell time-steps the more effective is the temporal adaptation; as is true for 
spatial adaptation. 

The finite volume, node based, scheme permits independent integration 
for each cell but requires increased artificial dissipation for finer cells at the 
spatial interfaces when multiple grids are used for steady state situations. 
The scheme also allows for separate pointer systems for spatial and temporal 
adaptation procedures. 

The spatial data base tracks the spatial level, parent cell, and the sur- 
rounding nodes of each cell in the domain. Similarly, information about 
cells adjacent to each node must be known as well as that for special nodes 
such as the boundary and spatial interfaces. The temporal data base tracks 
the number of cells and the sequence of integration during each time-stride. 
This pointer system must be updated after each time-stride for assignments 
of time-steps, determination of the cell temporal level and its allocation into 
a cluster at the same level, determination of nodits, and constraining of time- 
steps across them (allowing only factors of two). Some of this represents an 
overhead but when compared to the gain achieved in efficiency proves to be 
well worth doing. Depending upon the problem, the spatial data base updat- 
ing may not be required as frequently as the temporal data base. For steady 
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state problems the number of changes in the spatial pointer system equal 
the number of spatially embedded levels desired and the adaptation can be 
performed at either specified iteration intervals or at any specified residual 
level. Similarly, for unsteady problems in which the characteristic feature 
speeds are relatively small the adjustments to the spatial pointer system 
are infrequent. However, when high feature speeds arise, either the time- 
stride size must be kept small or the spatially embedded clusters enlarged, 
so that the features do not move out of their respective clusters during a 
given time-stride. Furthermore the actual cell time-steps may change during 
a time-stride and hence the nodits might have to be moved to regions involv- 
ing relatively less gradients. The process of enlarging spatially embedded 
clusters and translation of nodits can become computationally expensive; a 
balance is required between the competing effects of temporal level stiffness 
and temporal adaptation. The enlarging of the spatially embedded grids 
is accomplished by first determining the contiguous clusters of the refined 
cells and then refining those coarse cells which are outside and adjacent to 
be identical in level to those just inside the boundary, and repeating this 
a specified number of times. The translation of nodits similarly determines 
the nodits of clusters for contiguous cells with the same temporal level, and 
then halves the time-steps of the neighbouring cells with the larger time-step, 
repeating a specified number of times. 

For all sample cases the numerical solutions with adaptation procedure 
were comparable in accuracy to globally fine grid solutions, and are in good 
agreement with previous works. Since the savings in CPU time are sub- 
stantial for this quasi-one-dimensional study, it does appear promising to 
introduce temporal adaptation concurrently with spatial adaptation for mul- 
tidimensional, unsteady reacting flow fields. 
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